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We consider the possible mechanical instabihty of an ultracold Fermi gas due to the attractive 
interactions between fermions of different species. We investigate how the instability, predicted by 
a mean field calculation, is modified when the gas is trapped in a harmonic potential and quantum 
effects are included. 
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I. INTRODUCTION 

The field of ultracold Fermi gases has recently seen a very strong development. After reaching the degenerate 
regime, experiments have dealt with mixtures of two hyperfine levels [1-3] . This has been done mostly in the vicinity 
of a Feshbach resonance, which allows quite conveniently to adjust the effective interaction. Indeed by changing the 
magnetic field, one can start with a small negative scattering length, make it more negative, let it have a jump from 
— oo to +CX3 at the resonance, and then have it decrease to small positive values. A major purpose of these experiments 
is to look for the BCS transition, which should occur in particular for negative scattering length a. Naturally one 
expects the critical temperature to be higher for larger |a|, since this corresponds to a stronger attractive interaction. 

On the other hand one may expect that this overall attractive interaction may give rise also to a collapse instability. 
This would be similar to the one very much studied in Bose-Einstein condensates [4] , where this instability prohibits 
the formation of condensates with a large number of atoms. It is reasonable to think that similarly the BCS instability 
would be in competition with a collapse. The first work where the BCS transition for ultracold atoms was explicitely 
considered [5] studied also this collapse within the mean field approximation, considering in particular the dependence 
on the ratio between the number of atoms in the two atomic population expected to form Cooper pairs. The vicinity 
of this collapse would also be particularly favorable [6] to the BCS transition. 

More specifically let us consider only the case of two equal populations of atoms corresponding to two different 
hyperfine levels, and restrict ourselves to the T = case. A scattering length a corresponds to an effective interaction 
constant g = Anh^a/m between unlike atoms, giving rise to a mean field contribution gn/2 to the chemical potential 
where n is the total atomic density. Including kinetic energy the total chemical potential is given in this 
approximation by fJ-in) — h'^kp/2m — \g\n/2 with Stt^u — kp. The instability is obtained from dfi/dn — which gives 
for the critical density A = 2kp\a\/Tr — 1. 

Quite remarkably recent experiments [1-3,7] on mixtures of fermionic atoms with two different hyperfine states 
have not observed this instability when they have been through the Feshbach resonance, where the scattering length 
becomes infinitely large. Since these experiments have reached quite low temperatures it is very unlikely that they 
missed this transition because they did not go at low enough T. On one hand the absence of this instability is quite 
fortunate since it allows experiments to reach all the possible range of scattering lengths, without any limitation. In 
particular this gives the possibility to produce molecules, which has been extremely fruitful very recently. Nevertheless 
the failure of the mean field calculation [5] is somewhat striking since, although there is no reason to believe mean field 
to be quantitatively correct, it gives quite often reasonable qualitative estimates. Note however that there is nothing 
basically wrong in the idea that a collapse instability should exist since it has been indeed observed experimentally 
[8] in a fermion-boson mixture with an effective attraction between fermions and bosons. 

In this paper we investigate a possible source for this disagreement between this simple mean field estimate [9] 
and experiments. An indication in this direction can be found in the results of an investigation of the hydrodynamic 
modes in harmonic trap [10] . The hydrodynamic framework implies that the quantum length scale does not appear 
in the study, except indirectly in the equation of state of the dense gas. The surprising result is that no mode with 
zero frequency is actually found, even when the system reaches the density for collapse at the center of the trap. All 
the modes are found to have non zero frequencies, whereas one would expect the collapse instability to manifest itself 
by the appearance of a mode with zero frequency, as it is found for a Bose condensate. Indeed one can show, from the 
starting equations, that it is a general feature of these hydrodynamic modes to have non-zero frequencies [11]. This 
result may suggest that, in a trap, the collapse may be missed in some way. 
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In order to explore this problem and go beyond the macroscopic scale by including quantum effects, we perform in 
this paper a semiclassical microscopic calculation to find if there is a zero frequency mode. Since experiments have 
already reached very low temperatures, and that this is the most favorable situation for appearance of the instability, 
we restrict our exploration to the T = case. We find indeed that there is one. The lack of zero frequency mode in 
the hydrodynamic framework has actually a simple physical explanation. The gas becomes unstable when the density 
at the center of the trap reaches the critical density. However the gas in all the other parts of the trap has a density 
lower than this critical density. So the overall system is, so to speak, not soft enough to a have a zero frequency mode. 
When one goes to a microscopic calculation, the unstable region, which had at the macroscopic scale a zero extension, 
gets a finite extension of order of the microscopic scale. In this way this region can produce a strong enough softening 
and give rise to a zero frequency mode, which we obtain below explicitely. 

In the next section we begin our investigation by making a simple RPA approximation to find the microscopic effect 
of interactions. Our treatment takes into account the modification of the density distribution due to the interactions. 
This allows us to find out explicitely the important features in this problem, and in particular to point out the 
relevant length scales. We can then generalize our approach by getting rid of the RPA approximation, and show how 
to deal with the problem on the quite general grounds of Fermi liquid theory. We find indeed that the threshold for 
instability is modified by the trap and obtain explicitely this modification, as well as the shape of the mode responsible 
for the instability. However, although this modification could be sizeable for the very anisotropic traps used in some 
experiments, it appears unlikely that it is responsible for the overall disappearance of the instability, as it is observed 
experimentally. 



II. THEORETICAL TREATMENT OF THE INSTABILITY IN A TRAP 



We consider now the above atomic gas of fermions with mass m, in an isotropic harmonic trap of frequency fi, 
giving rise to the harmonic potential V{x) = ^mfPr"^ with = x^. An (undamped) eigenmode corresponds to an 
infinite response of the system excited at the frequency of the mode. Since we are interested in a zero frequency mode 
[12], we have to consider in addition a static perturbation 6V{x), which will induce a static density fluctuation (5p(x). 
The collapse instability will correspond to a divergent density fluctuation. The linear response theory gives [13] : 

5p(x) dx'H(x,x')(51^(x') (1) 

where H(x, x') = —j-f dt6{t){[p{x,t), p{'x' ,0)]) is the zero frequency density-density (retarded) response function. 
Therefore an instability corresponds to a divergent eigenvalue of the kernel operator H(x, x'), which is real symmetric. 
For interacting particles, one has to make use of some kind of approximation to calculate H. The simplest one, which 
will reduce as we will see to the mean field result for an infinitely wide trap (in which case the system would be 
homogeneous) is the Random Phase Approximation (RPA) and this is the one we will use here. As it is well known, 
it is equivalent to sum up an infinite series of bubble diagrams. Since the fermions we deal with interact with a 
very short range atomic size potential, we can use a contact potential U{x) = gS{x) , with g — 'inh^a/m, for the 
interaction potential, where g < since we consider an attractive interaction. In this case the RPA for our trapped 
fermions reads: 



H(x,x') = H°(x,x') +5 j dxi H°(x,xi)H(xi,x') 



(2) 



This Eq. (2) reads formally H = H° + gH°H, and its formal solution is H = (1 - .gn°)"^H". This shows that the 
eigenvectors of H and H° are the same, and that the instability we look for appears when the smallest (negative) 
eigenvalue of H° equals 1/g. The general eigenvalue equation for H° (eigenvalue a, eigenvector ip) reads: 



j dx' H°(x, x') ^(x') = a (^(x) 



(3) 



It is then convenient to introduce the Wigner transform H^(q, R) = J dre "i ''H°(R + r/2, R — r/2). Eq.(3) then 
becomes: 

Then we use the fact that there is a large number of trapped particles or equivalently that the chemical potential 
is much larger than the level spacing %Q, . This allows to make use of a semiclassical treatment by considering that 
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the trapping potential is slowly varying. In this case we can use for n^(q, R) its homogeneous value, evaluated with 
the local value of the particle density. We therefore make the approximation n^(q, R) f« n°(q), where 11° (q) is the 
response function of the homogeneous system with a Fermi wave vector kpiR)- The local Fermi wave vector is related 
to the equilibrium density of the cloud n{R) = kp/Sw^, determined by the equation: 

M(n) + 1/2 mO^r^ =/i (5) 

where /u(n) is the chemical potential and jl is the overall chemical potential. 

Our next step is to take advantage of the length scale d we expect physically for the instability mode we are 
interested in. Clearly the instability will occur at the center of the trap since this is where the local particle density 
is highest. On the other hand this mode is a collective phenomenon involving a large number of particles, so it must 
occur over a typical scale large compared to the interparticle distance, which is itself of order kp^. This leads us to 
look for a mode which satisfies fc^^ <C d <C i?o- This relation implies that the typical wavevectors entering the Fourier 
expansion of </9(x) are small compared to kp- From Eq.(4) it is then seen that the wavevector q in n^(q, R) must also 
be small compared to kp- This allows us to expand 11° in powers of q. Since n°(q) is just the free particle response 
function, we have [13] expHcitely n"(q) w ~^^mkF{l — j2Q^/kp)- When we insert this expression in Eq. (4) and 
perform the integrals, we find the following second order partial differential equation for the density fluctuation tp 
corresponding to the instability mode: 

Aip + krVi^^.V^fi + \kF^pA{^) + 12(fc| + 27r^— a)ip = (6) 
kp 4: kp m 

If we consider now the order of magnitude of the three first terms of Eq.(6), we notice that the second and the 
third term contain derivatives of kp, while the first one contains only derivatives of Since the length scale Rq for 
the variations oi kp is much larger than the length scale d for the variations of the second and third terms are 
negligible compared to the first one. 

In order to solve explicitely Eq.(6) we consider as a first step, in the following subsection, the simple case where 
the modification of the density distribution due to the interactions is not taken into account. This will allow us to see 
clearly the relevant length scales. This simple case is equivalent to assume the free particle relation /^(n) = fi^ kp/2m 
for the equation of state. We will then take consistently into account interactions in kF{R) in the next subsection. 



A. Simple case 

In this simple case we have merely [14]: 

fcf(i?) = 4(i-i?Vi?2)i/2 (7) 

where the Thomas- Fermi cloud radius Rq and the maximum Fermi wave vector kp are related by [14] {k%f /2m = 
(1/2) mfi^ -Rq. Both are directly related to the particle number N in the trap. Coming back to Eq. (6), since again 
d <C Rq, kp{x) is a slowly varying function in the considered domain of a; = ]]x]] ~ d. We can therefore expand kp{x) 
up to second order around x = 0, using Eq.(7). Setting the eigenvalue a = 1/g at the instability, we can introduce 
the coupling constant A = —^^fnkpg = —^kpU, and we find: 

+ 12 - l/A) - {k°p/Rof{l - l/2X)x^] ip = (8) 

Going up to fourth order in the expansion of kpix) would give a term of order (fc^)^/(i?o)^ a;^, which is smaller than 
the term in Eq.(8) by a factor (x/i?o)^ ^ 1- It is therefore justified to stop the small x expansion to second order 
to get Eq.(8) from Eq.(6). 

Now Eq.(8) is the Schrodinger equation for the 3D harmonic oscillator of frequency w for a state of energy E, 
provided we set ?i = m = 1, together with: 

£; = 6(4)^(1-1) (9) 

and 

= 12(4/i?o)'(l - h (10) 
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Since the instability mode we are looking for is naturally localized in the center of the trap, we arc looking for the 
bound states of this harmonic oscillator. The critical value of the coupling constant A is directly obtained from 
the energy E of the oscillator by Eq.(9). Since the first instability will occur for the smallest value of the coupling 
constant, we are looking from this equation for the smallest possible value of the energy. We check that we recover 
properly the homogeneous case by taking the limit Rq oo. Indeed in this case the harmonic oscillator frequency lo 
goes to zero, which implies that all the bound states energies go to zero. From Eq.(9) this gives A = 1 as expected. 
Coming back to the trapped case we will find the lower energy among the isotropic s— wave states. We have for 
these states the quantization relation E = w(2n + 3/2), with n = 0, 1, • • •. For a given n, this determines A at the 
instability. The smallest value of A is obtained for n = 0, corresponding to the gaussian mode (p{x) = e~2"^ . This 
yields 1 - 1/A = a/3/4(1 - l/2A)(i?ofc^)"^ < 1. This result shows that we have A ~ 1. Therefore to first order in 
1/Rokp we have at the instability: 

Naturally we have to check the consistency of our calculation by looking at the size d of the instability mode. Since 
we have A ~ 1, we can just set A = 1 in the equation Eq.(lO) for w^. We get u = \/6kp/Ro leading to a gaussian mode 
(fix) = exp{~^x'^/(P) of width d = = 6-^/^{Ro/k"py^^ = Except for the numerical coefiicient 

this is the size of the single particle ground state in the harmonic trapping potential. This result for the width of the 
mode is completely consistent with our starting hypotheses fc^^ <C d <C -Ro- 



B. Self-consistent calculation 



Now, since we have the quantitative situation under control, we come back to a consistent description of the density 
distribution in the atomic cloud, taking interactions into account in the calculation of the Fermi wave vector. Taking 
for the chemical potential the Hartree approximation, we have — fi^kl^ /2m — \g\n/2 with n = kp/Sn^. When 
A = 1, the static compressibility {ndii/dn)~^ diverges at the center of the trap. As a result when we take the derivative 
of the chemical potential with respect to the Fermi wavevector, we find that it is zero. This implies that kpiR) is a 
linear function of R close to the trap center, instead of being quadratic as in Eq.(7). As we have seen in the preceding 
subsection, we need to know kpiR) only close to the center. This is easily done and one finds: 

k% ~' V3Ro' ^ ' 

where Rq is again the Thomas-Fermi radius of the atomic cloud. But it is now related to the Fermi wavevector at the 
center kp by hkp = VSmflRo, instead of hkp = mflRo as above. Following the same procedure as before, we insert 
this expression for the density into Eq.(6). Keeping only the dominant terms, we get: 

A^ + 12(4)2(l-l-i=-^V = (13) 

For s-wave solutions this equation can be reduced to Airy function differential equation: 

i^"{y) = {y-yo)^{y) (i4) 

provided that we rescale the position x according to a; = Dy, with the new length scale D given by D/Rq = 
{4y^{k°pRof))-^/^, and we introduce the new function ip = X(p. We find yo = 12 {k°p Df (1 - 1/A). We note that 
the power law dependence of D/Rq on kpRo is slightly different from the one we have found for d/Ro. Nevertheless 
our starting hypotheses kp^ <C d <C i?o are still satisfied. The boundary conditions V-'(O) = and ^"(+00) = impose 
that 'ip{y) = Ai{y — yo) where yo » 2.3 is the first zero of the Airy function. We finally get for the coupling constant 
at which the instability arises: 

This result is similar to the one we have found above Eq.(ll) for our simple case. However the dependence on kpRo 
is somewhat weaker since the exponent is 2/3 instead of 1. 
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Our result Eq.(15) is coherent with what might be expected physically. Indeed this can not be the density right 
at the center of the trap which is relevant for the instability. One has rather to consider the average density over a 
region of typical size a few kp^. This average density is lower than the nominal density right at the center. We expect 
therefore that the threshold for the instability is raised, compared to what one could obtain by considering only the 
density at the center. This is just what we obtain. Nevertheless the shift of the instability is rather small since it 
is of order l/{Rokp)^/^ . On the other hand this is coherent with the fact in the limit Rq ^ oo we have to recover 
the homogeneous case with no shift at all. Hence the result has to depend on the ratio between the microscopic 
length kp^ and the macroscopic one Rq, wich is small. Therefore we come to the conclusion that the effect of the 
trapping potential is unlikely to be responsible for the lack of collapse found in experiments. On the other hand it is 
not completely clear that the very elongated shape of the traps used in most of this kind of experiments would not 
play some role. Indeed in this case the above ratio would not be that small. Nevertheless we have not explored this 
more complex situation. However in at least the ENS experiment [15] the trap is not so far from isotropic, so strong 
anisotropy does not hold for explaining the lack of collapse. 



C. General case 



Finally we generalize our treatment by getting rid of our approximate evaluation of the static response function of 
the system. Indeed we used above the simple RPA to obtain it. However our above derivation makes it clear that, 
because of our scmiclassical approximation, all what we need is the response function of the homogeneous system. 
Hence we can use the general framework of Fermi liquid theory to discuss it. Let us recall that this framework is an 
exact one, and that it applies to our case of interest, that is a strongly interacting neutral Fermi system [16], just as it 
does for normal and superfluid liquid ^He. Naturally there is a counterpart to the fact that this framework is exact, 
which is that it does not provide the explicit values of the constants it introduces. However this is not so important 
here since our approximate treatment provides us already with order of magnitude for these constants. On the other 
hand it is quite important to know that the only things which are not exact in our treatment are the values of these 
constants. 

We found above from the RPA that the response function H for the interacting system could be obtained from the 
response function n° of the non interacting system by H~^ = (n°)~^ — g. Now in Fermi liquid theory [16], for zero 
wavevector, the exact value of H is given by: 

where A^o = m*kF/TT'^, is the density of states, with m* the effective mass, and Fq the £ = symmetric Landau 
parameter. We see that we have merely to replace the bare mass by the effective mass, and replace g by Fq/Nq to 
obtain the exact result (note that we have a factor 2 difference with our above explicit expression for the density 
of states because we consider the two spin populations). In particular the collapse instability corresponds to the 
well-known condition Fq = —1. Hence we only need to know the exact dependence of A^o and Fq as a function of 
particle density in order to recast our above calculation under an exact form. More precisely we need only to know 
this dependence in the vicinity of the collapse density. Finally in order to write the generalization of Eq.(13) we need 
not only our response function for zero wavevector, but also the lowest correction due to the fact that we work at 
nonzero wavevector q. For dimensional reasons this correction will be in (q/kp)^ but the coefficient will not be given 
by the free gas result, as we have done above. More generally this coefficient is beyond the reach of standard Landau's 
Fermi liquid theory. This leads us to write finally the small q expansion; 

where we had b = 1/12 in the free particle case. Hence the framework we used in our calculation is exact, only the 
constants which come in have to be modified. We do not rewrite here the generalization of our above calculations 
since it is straighforward to do it and the results are not expected to lead to qualitative changes with respect to our 
above ones. 



III. CONCLUSION 



In this work, we have considered the zero frequency unstable density fluctuations for ultracold fermions in a 
harmonic trap. More precisely we have studied how quantum effects modify the corresponding mode. We have used 
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a semiclassical treatment justified by the fact that the trap is large compared to the microscopic quantum scale. In 
a first step we have treated the interactions within the RPA approximation. Wc have then made a fully general 
analysis within Fermi liquid theory. Our results show that there is a zero frequency mode, with a size large compared 
to the inverse Fermi wave vector at the center of the trap but small compared to the cloud size. The threshold for 
instability is modified by the trap. However, even though we find a modification in the threshold instability, this 
docs not seem to explain the absence of instability, as it is observed experimentally. It is interesting to note that, by 
contrast, quantum effects have a strong infiuence on the corresponding instability for Bose systems. This is basically 
because in this case the correspondent of our dimensionless parameter Rokp becomes of order unity. It would be of 
interest to devise a model allowing to go continuously from the Fermi case to the Bose case, and see the importance 
of quantum effects grow along the way. 

* Laboratoire associe au Centre National de la Recherche Scientifique et aux Universites Paris 6 et Paris 7. 
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